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ABSTRACT 


Parabolic  equation  models  solved  using  the  split-step  Fourier  (SSF)  algorithm,  such  as 
the  Monterey  Miami  Parabolic  Equation  model,  are  commonly  used  to  predict 
underwater  sound  propagation  in  deep  and  shallow  water  environments.  Previous  studies 
have  shown  that  the  SSF  algorithm  is  very  accurate  in  shallow  water  when  there  is  no 
density  discontinuity  between  the  water  column  and  the  sediment,  but  less  effective  in  the 
presence  of  realistic  density  discontinuities  due  to  phase  errors  that  accumulate  after  a 
few  kilometers.  In  this  thesis,  the  standard  density-smoothing  approach  and  an  alternative 
hybrid  split-step/finite  difference  method  are  compared.  The  goal  is  to  decrease  the  phase 
errors  and  increase  the  model’s  long-range  accuracy.  Different  depth  meshes  and  range 
step  sizes  are  implemented  in  the  algorithm  to  find  the  optimum  results  for  both 
approaches.  It  is  shown  that  the  density-smoothing  method  provides  better  results  with 
small  range  step  sizes,  while  the  hybrid  method  is  more  effective  using  longer  range  step 
sizes.  However,  the  smoothing  approach  provides  a  more  stable  convergence  of  results, 
whereas  the  hybrid  method  solution  is  more  sensitive  to  change  in  computational  grid 
sizes.  A  more  detailed  examination  of  the  density  smoothing  approach  suggests  good 
accuracy  for  a  few  kilometers,  while  the  hybrid  method  provides  improved  agreement 
with  a  benchmark  solution  at  longer  ranges. 
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I.  INTRODUCTION 


A.  BACKGROUND 

The  study  of  the  parabolic  equation  (PE)  approximation  to  the  elliptic  Helmholtz 
wave  equation  dates  back  to  mid- 1940s,  when  Leontovich  and  Fock  introduced  the  PE 
method  to  the  problem  of  radio-wave  propagation  in  the  atmosphere.  Since  then, 
propagation  modeling  using  the  PE  has  been  successfully  applied  to  microwave 
waveguides,  optics,  laser-beam  propagation,  plasma  physics,  seismic  propagation  and 
underwater  acoustics  [1].  Tappert  was  the  first  to  introduce  the  PE  method  for  underwater 
acoustic  propagation  in  1974  [2],  and  the  split-step  Fourier  (SSF)  algorithm  to  solve  the 
PE  was  developed  by  Hardin  and  Tappert  in  1973  [3],  Since  then,  many  articles  have 
been  published  to  improve  PE  propagation  modeling  in  the  ocean.  According  to  a  survey 
made  in  1990,  over  a  period  of  15  years,  more  than  120  articles  and  technical  reports 
related  to  PE  developments  in  underwater  acoustics  had  been  published,  as  well  as  many 
more  after  that  time  [1]. 

The  main  advantage  of  the  parabolic  equation  is  modeling  underwater  sound 
propagation  in  the  ocean  at  long-range  and  in  range-dependent  environments  [4].  Several 
methods  and  algorithms  have  been  introduced  to  solve  the  parabolic  equation.  The  most 
common  ones  are  the  SSF  algorithm,  the  implicit  finite  difference  (IFD)  algorithm,  and 
the  finite  element  (FE)  method.  There  are  many  advantages  and  disadvantages  to  these 
techniques.  For  SSF,  the  primary  advantage  is  its  speed  and  simplicity  [5].  Also,  the  SSF 
is  efficient  for  long-range,  narrow-angle  propagation  problems  in  range-dependent  media. 
In  environments  causing  strong  bottom  interactions,  however,  it  loses  some  of  its 
accuracy.  IFD  and  FE  methods  give  more  accurate  results  in  wide-angle,  bottom¬ 
interacting,  and  range-independent  media.  However,  they  are  generally  less  efficient  and 
less  stable  than  the  SSF  algorithm  [1]. 

The  standard  parabolic  approximation  originally  developed  by  Tappert  was  valid 
only  for  small  angles,  and  its  accuracy  was  degraded  by  the  phase  errors.  Several 
attempts  were  made  to  reduce  these  errors  and  support  wider  angles  of  propagation. 
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Thompson  and  Chapman  introduced  the  wide-angle  parabolic  equation  (WAPE) 
approximation,  which  was  based  on  Feit  and  Fleck’s  operator  splitting  [6].  This  approach 
treats  the  square -root  operator  in  a  way  that  differs  from  the  standard  parabolic  equation 
(SPE).  It  improved  the  parabolic  equation  by  increasing  the  propagation  angle  from 
around  15  degrees  to  40  degrees,  reducing  the  phase  errors  while  retaining  the  usage  of 
the  efficient  SSF  algorithm  [4], 

As  previously  stated,  the  SSF  algorithm  is  known  to  be  less  accurate  in  the 
presence  of  strong  bottom  interactions.  In  order  to  treat  the  density  contrasts  at  the 
interfaces  between  differing  media  (e.g.,  water  and  sediment),  a  general  smoothing 
function  is  applied  to  the  interface.  This  smoothing  approach  is  the  primary  cause  of  the 
phase  errors  due  to  density  discontinuity  at  the  interface  for  long-range  propagation.  In 
1996,  Yevick  and  Thomson  introduced  a  hybrid  split-step/finite  difference  (SSF/FD) 
technique  for  modeling  acoustic  propagation  in  the  presence  of  non-uniform  density 
profdes  [7].  This  method  offers  a  practical  way  to  adapt  existing  SSF-based  models  by 
separating  the  operator  approximations  into  density-dependent  and  density-independent 
components  rather  than  relying  on  smoothing  parameters.  An  overview  of  this  method  is 
presented  in  this  thesis  with  specific  numerical  examples  showing  the  improvements  in 
the  solutions. 

B.  PROBLEM  STATEMENT 

The  Monterey-Miami  Parabolic  Equation  (MMPE)  model  was  developed  in  the 
mid-1990s  and  since  then  has  been  tested  for  several  existing  benchmark  scenarios.  The 
MMPE  method  computes  the  underwater  acoustic  propagation  by  using  the  SSF 
algorithm  [8].  There  are  various  versions  of  MMPE,  including  both  2D  and  3D  versions, 
versions  that  incorporate  rough  surfaces,  oceanographic  internal  waves,  and  a  variety  of 
other  environmental  perturbations.  In  all  cases,  however,  it  applies  the  WAPE 
approximation  of  Thomson  and  Chapman  [4]. 

Various  techniques  have  been  developed  to  treat  the  density  discontinuity  in  the 
bottom  interaction  in  order  to  increase  the  accuracy  of  MMPE.  In  this  thesis,  our  goal  is 
to  compare  the  effects  of  Tappert’s  standard  density  smoothing  approach  and  Yevick  and 
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Thomson’s  hybrid  split-step/finite  difference  algorithm  [7]  in  the  case  of  the  bottom 
interactions  for  shallow  water.  We  will  perform  the  propagation  in  a  simple  Pekeris 
waveguide  defined  by  Paul  Hursky  [9]  and  previously  investigated  by  Owens  [8]  and 
Smith  et  al.  [10]. 

In  order  to  avoid  any  issues  with  implementation  details  in  the  existing  MMPE 
model,  a  Matlab  version  of  the  SSF  algorithm  has  been  developed  for  this  thesis  based  on 
the  same  operator  approximations  as  the  MMPE  model.  It  was  also  written  only  for  the 
simple,  range-independent  Pekeris  waveguide  that  motivated  this  study. 
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II.  THEORY 


A. 


PARABOLIC  EQUATION  DERIVATION 


We  start  with  the  2-D  Helmholtz  wave  equation  in  cylindrical  coordinates, 


]_8_ 

r  dr 


dp 

dr 


d 

+  p — 

dz 


1  dp 
p  dz 


+  kln  p  =  0, 


(1) 


where  kn  =  —  is  the  reference  wave  number,  n  =  — 5 —  is  the  acoustic  index  of 
c0  c( r,z) 

refraction,  c0  is  the  reference  sound  speed,  c(r,  z)  is  the  spatially  varying  acoustic  sound 

speed,  and  p(z)  is  the  density  that  is  assumed  to  change  only  with  depth.  All  features  of 

the  environment  are  presented  in  c(r,  z)  and  p(z) .  We  define  the  PE  field  function  VP 

according  to 


'P (r,  z)  =  y[k~rp(r,  z)  e  'v , 
which  leads  to  the  fundamental  parabolic  equation 

PVT/ 

dr 

where  the  square -root  operator  Q  is  defined  as 


„  old 
Q  =  iln~  +  ^-p — 

V  kr  dz 


' i_y 

yP  dz  y 


We  can  also  define  the  field  function  as 

which  satisfies  the  PE 


(2) 


(3) 


(4) 


(5) 


dr 


=  ik0(Q- 1)*. 


(6) 
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Here,  the  square-root  operator  takes  the  more  compact  form  of 


0  = 


i  a2 

kl  dz2 


(7) 


and  the  effective  index  of  refraction,  n  ,  is  defined  as 


~2  2  , 
n  =  n  +  - 


2  kl 


i  d2P 

p  dz2 


1  dp 2 
\P  dz  , 


(8) 


The  original  SSF  algorithm  is  based  on  (8),  which  is  discussed  further  in  the 
section  describing  the  smoothing  approach.  Either  PE  allows  for  a  one-way  marching 
algorithm  solution  of  the  form 

T(r+Ar,z)  =  exp[/k0Ar(g-l)]fi/(r,z)  (9) 


or  similar  for  'P  and  Q .  We  can  define  the  operators  as 


2  ,  132  , 

s  =  n  -1,  p  =  ——,  and  y  ■■ 
^  dz 


The  original  square  root-operator  then  becomes 


1  dp  d 

kip  dz  dz 


2  1  d 


Q  =  <  ir+PiP— 


f  l  d  A 


kl  dz 


p  dz 


^  1  dp  d  1  5’  ^ 

p1  dz  dz  p  dz2 


=  ]j  l  +  (n2-l)  + 
=  yjl  +  s  +  p  +  y 


15“  1  dp  d 


kl  dz~  k0p  dz  dz 


(10) 


(11) 


which  requires  further  approximation  in  order  to  implement  the  marching  algorithm. 

The  MMPE  model  utilizes  the  Thomson-Chapman  WAPE,  defined  by  the 
operator  splitting  of  Feit  and  Fleck  [6], 

g;  =  >/lT7  +  A/lT^  +  >/rKF-2.  (12) 
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Then  the  propagator  in  the  marching  algorithm  takes  the  form 


exp  [ik0Ar(Q'  - 1)]  =  exp 

iS ^-\/l  y  H"  yj  1  +  fU  +  -\J  1  -\-  £  —  3  j 

«  exp 

iS^l  +  y-l^j 

exp 

iS^l  +  ju  -l) 

exp^/A^Vl  +  £  - 1  j 

=  exp 

is[4^+r- i) 

exp 

exp[z(7(n-l)] 

(13) 


where  S=koAr  and  the  non-commuting  terms  of  the  operators  e,  jx,  and  y  are  neglected  in 
the  second  line.  As  strictly  noted  by  Yevick  and  Thomson  [7],  the  operation  order  in  (13) 
is  very  important,  and  the  y  term  should  be  applied  at  the  end  of  the  range  step, 
following  application  of  the  n  tenn. 

The  operator  currently  used  in  MMPE  does  not  contain  a  y  term,  since  the 
density  tenns  are  built  into  the  effective  index  of  refraction.  In  this  case,  we  simply 
define 

a=Vrr^+>/m-i,  (i4) 


where  £  =  n2 -1.  The  WAPE  approximation  does  not  consider  the  lowest  order  cross 
tenns  between  operators,  so  the  lowest  order  correction  becomes 


a=a- 

=  Q[- 


]_ 

8 

]_ 

8 


-  (s  +  /j  +  y)  - s1  - /u 2 

[^(//  +  y)  +  (//  +  y)^  + 


w  +  w\ 


(15) 


For  completeness,  note  that  the  lowest  order  correction  to  the  smoothing  approach  would 
be 


Qi=Qx- 

=  6- 


]_ 

8 

1_ 

8 


( s  +  fif-s 2 
[£•//  +  fxs\. 


(16) 


In  the  SSF  algorithm,  a  and  u  are  just  scalar  operators  in  the  physical  (z)  and 
wavenumber  spaces  ( kz ),  respectively. 
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If  the  y  terms  and  cross  tenns  are  added  later  as  corrections,  then  we  can  define  an 
intennediate  step  solution  as 


VF  (r+  A  r,  z)  =  FFT 


' 

f  1 - TT  > 

iS 

1  k- 

exp 

J  ko  y 

_ 

_ 

■  FFT  [exp  (iS(n- 1))' ^(r,  z)] 


,  (17) 


where  we  have  replaced  y]l  +  ju=>  1 — f  and  y/l  +  s  =>n,  consistent  with  the 

v  K 

application  in  the  ^-domain  and  z-domain.  Note  that  the  complete  range  step  in  the 
existing  MMPE  code  is  obtained  by  replacing  n  =>  ft ,  where  T(r,z)  is  replaced  by 

vF(r,z). 


B.  DENSITY-SMOOTHING  APPROACH 


The  treatment  of  the  density  discontinuity  in  the  MMPE  model  utilizes  a  density¬ 
smoothing  approach  at  the  water-bottom  interface.  As  presented  in  (8),  instead  of  the 
standard  index  of  refraction,  n,  we  define  an  effective  index  of  refraction,  h ,  by 


~2  2  , 
n  =  n  +  - 


2  kl 


1  d2p 

p  dz2 


1  dp2 

\P  dz  j 


(18) 


The  density  discontinuity  at  the  water-bottom  interface  can  be  defined  as 

P(z)  =  Pw  +  (Pb  ~PJ  H(z-Zj), 


(19) 


where  pw  and  pb  are  the  water  column  and  bottom  sediment  densities,  respectively,  and 
are  assumed  to  be  constants.  The  function  FI  serves  as  an  approximation  to  the  Heaviside 
step  function  and  defines  a  smooth  transition  according  to 


c 

1  +  tanh(-Z— ) 
2L 


(20) 


where  <f=  z-z/„  and  L  defines  the  mixing  length. 
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The  first  and  second  derivatives  of  the  density  are  defined  in  the  effective  index  of 
refraction  as 


dp_ 

dz 


dz  4  L 


seclr 


2  L 


(21) 


and 


d2p  .  ,  d2H 

—  =  (Ph~Pj - 


dz1 


(Pb-Pj 
4  L2 


tanh 


f 

—  sech2 
2  L 


C 

v2  L) 


(22) 


These  expressions  are  then  implemented  into  (18)  to  define  the  effective  index  of 
refraction. 


C.  HYBRID  SPLIT-STEP  /  FINITE  DIFFERENCE  IMPLEMENTATION 

TECHNIQUE 

The  hybrid  split-step/finite  difference  PE  algorithm  for  variable  density  media 
was  introduced  by  Yevick  and  Thomson  in  1996  [7].  Instead  of  relying  on  a  smoothing 
function  to  transition  the  density  discontinuity  across  the  interface,  they  used  a  hybrid 
technique  to  specifically  solve  for  the  terms  containing  the  density  contrast.  This 
approach  can  easily  be  implemented  into  existing  SSF  codes.  To  accomplish  this,  they 
separated  the  effects  of  the  density  by  splitting  the  propagator  into  density-dependent  and 
-independent  components,  as  described  in  the  previous  section.  In  this  approach,  the 
density  propagator  is  solved  by  using  the  finite  difference  technique  and  a  Pade 
expansion  [7]. 


The  hybrid  approach  proposes  that  (Q[  requires  an  additional  operation  to  (17), 
given  by 


vF'(r+  Ar,z)  =  exp  id(^j\  +  y  -1  j 


Tfi+Ar,  z). 


(23) 


If  we  intend  to  invoke  only  the  Q[  approximation,  then  this  completes  the  range  step 
calculation,  and  T* '  =>  T .  Here,  the  additional  operation  is  approximated  by  a  Pade 
expansion  defined  as 


9 


1  +— (l  +  iS)y 

exp|  iSly/T+x-l)  * — y - . 

1+1(1-,% 


(24) 


In  order  to  employ  a  finite  difference  approach,  we  then  write  (23)  as 


'P  (r+  A  r,  z)  = 


l  +  i(H-0)y 


'P(r+  Ar,z) , 


(25) 


Since  direct  application  of  y  is  problematic,  Yevick  and  Thomson  proposed  an  alternative 
transverse  operator: 


p  d  (  1  d  ^ 
k\  dz\p  dz , 


(26) 


Then  y  =  //-//,  where  both  //'  and  //  are  well  defined  operators  on  Tf 


1.  Finite  Difference  Technique 

Consider  a  discretely  sampled  environment  and  associated  field  where  index  j=nb 
coincides  with  the  water-bottom  interface,  j<nb  is  in  the  water  at  the  shallower  depths 
and  j>nb  is  in  the  bottom  at  deeper  depths. 


Pi. )  <  nb)  =  pw; 

pQ  =  nb)  =  ^(pw  +  pb); 


(27) 


pQ  >  nb)  =  ph . 

p  d 

The  issue  here  is  the  evaluation  of  the  term  - 

kl  dz 


i  ay 

p  dz 


Assuming  a  smooth,  continuous  function,  a  3 -point  centered  difference  scheme  to 
evaluate  the  derivative  is  defined  by 


2  h 


(28) 
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where  f±l  is  the  function  evaluated  at  the  index  +1,  separated  at  a  distance  2 h,  to  give  an 
approximation  to  the  first  derivative  of  the  function  at  the  midpoint  index  0.  Equivalently, 


if  we  sampled  at  half  index  values,  then  /J  - 


1 

Ti 


When  the  term  f'  above  is  discretized  and  evaluated  with  this  approximation,  we 

obtain 

p  d(\  p _ i  _ 

k2  dzyp  dz  J  k\tSz  px  \  dz  )i  p  x  {  dz  )_i 

—  2  — —  2 

L  2  2 

— Of.-ToE— (f 

Pi  Pj 

2  2 

In  a  previous  analysis  by  Smith  [11],  implemented  by  Owens  in  2014  [8],  it  was  assumed 

that  since  p  is  discontinuous  at  the  interface,  then  the  factors  P-~^  c°uld  simply  be 

replaced  by  their  appropriate  values  above  or  below  the  index  being  evaluated.  However, 
it  has  been  detennined  that  for  the  finite  difference  approximation  to  hold,  p  must  be 
treated  as  a  continuous  function  in  the  vicinity  of  index  j.  In  this  case,  as  was  defined  by 
Thomson  [12], 


A  =\(Pi+P,),P  i  =-(p0+p_l). 
2  ^  2  ^ 


2.  Application  of  Finite  Difference  Approximations  to  Field  Operators 

Returning  to  our  field  operators  p  and  //' ,  we  find 


vr,  1  d2  1T(  1  d(8V\  1 

//T=  — ;rT=  — -  «— ; -  -  - 

kl  dz 2  kl  dz\dz  )  kl^z  \  dz  J\_  \  dz  )_}_ 

L  2  2  _ 


W-fot-CPo 


(31) 
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and 


kl  dz 


Po 


(  I  AT  ' 


Po 


\  p  dz  )  kl  (A  z)2 


Pi  P  , 


r  > 

1 

l  l 

1 

— ^1- 

—  +  — 

^0+ - ^  1 

Pl 

j  P,  P  , 

Pi 

2 

l  2  2j 

2 

(k0Az) 

Combining  (31)  and  (32)  withy  =  //'-// ,  we  obtain 

y  ¥  =  (//-//)¥ 


(32) 


where 


(£0Az)“ 


f  > 

f 

Po  y 

^1- 

Pi 

l  2  ) 

V 

A+  Al_2 

Pi  P  i 


\ 

f  > 

V  + 
T  o  T 

Po  j 
p  1 

V  1 

J 

l  2  J 

=77T?Ka-i)'i,,-(p-2)'i,»+(p--i)'I'-,] 

(£0Azj 


(33) 


P+-^  = 


Po 


+  A  I 
2  2 

Po 


p-  = 


(A+Po) 

Po 


2(A+P-i) 

P  =  P++P-  =  —  +  —  =  2p0 
Pi  Pi 

2  2 


1  1 

-  +  - 


(Pi+Po)  (Po+P-l) 


(34) 


Here,  the  definition  of  p  matches  the  definition  of  the  variable  p0  in  Yevick  and 
Thomson’s  paper.  However,  we  reserve  p0  to  designate  the  values  of  density  at  the 
central  index  j=0. 
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In  order  to  find  the  solution,  we  can  solve  the  tridiagonal  matrix  expression 


AY  =  BY , 


(35) 


where  A 


and  B 


i+i(i+i% 


.  A  and  B  are  just  identity  matrices 


except  in  the  vicinity  of  z  ~  zb,  where  y  is  not  zero.  Note  that  the  field  vectors  T  and  'F 

in  (25)  are  now  explicitly  written  as  vectors  over  depth,  T  and  T  .  The  solution  of  (35) 
completes  the  range  step  by  applying  the  density  operator  tenn  in  the  Qx  operator. 


3.  Higher  Order  Correction  Implementation 

Yevick  and  Thomson  [7]  also  proposed  a  higher  order  correction  term  in  Q'2 
given  by 


exp  -f  [e(n +/)+ (// + r)s + nr + m } | • 


(36) 


Since  //  +  /  =  //',  this  may  be  approximated  by 


f  jx  1  1  -  i '  iS(£ju'  +  ju's  +  /iy  +  y/u) 

exp  -T^'  +  ^  +  ^  +  W]  " - J - ■ 

'  l  +  —iS^£jU'  +  jU'£  +  jU/  +  /ju) 


(37) 


Note  that  if  y  =  0 ,  then  this  higher  order  correction  becomes  simply 


exp 


1-- ^—i5^s/u  +  /us) 
1  +  l^iS(£jU  +  /us) 


(38) 


This  tenn  can  be  applied  to  the  existing  MMPE  approach  based  on  the  operator 
approximation  Q2 . 


Yevick  and  Thomson  dropped  the  juy  and  ;///  terms  since  these  give  third-order 
derivatives  and  only  the  lowest  order  terms  are  used  in  the  finite  difference 
approximation.  Instead,  they  utilized  the  correction  as 
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exp 


^  1 4 - iS^sp'  +  p's) 


(39) 


This  produces  the  additional  correction 
'F(r+  Ar,  z)  =  exp 
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—  {sp'  +  p's) 


vF'(r+ Ar,z), 


(40) 


where  we  again  invoke  the  matrix  approach 


The  elements  of  the  matrix  A  and  B  are  defined  by  A  = 

.  Here,  splx¥  follows  (32)  as, 


!  +  ^(^'  + A) 


(41) 

and 


B  = 


+ A) 


£/PF  = 


^oPo 


(A:0Az) 


6 

— fiv 

1  1 

- + - 

Pl 

A 

P> 

2 

l  2 

2  7 

P  i 

2 


(42) 


(k0Az)‘ 


and  //(PT)  becomes 


Po 

8 

’  1  d(£¥)~ 

_  Po  5 

v 

k2 

K0 

dz 

p  dz  \ 

k02Az  <9z 

_Pv 

Po 


(k0Az)" 


Po 


(k0Az)~ 

1 

(k0Az)2 


V  2  2  2  27 


— 6.'!'.-^.)-— («0'P0-«_,'P_1) 
Pi  P  1 


—  ^1- 
Pl 


7  7 

1  1 

1 

—  + - 

r>  XT/  1  r-  VI/ 

^0  1  0  “  C-1  1  -1 

Pi  P  i 

Pi 

V  2  2  ) 

2 

(43) 


[p^i^i  -Po^o  +P-^-i'I/-i]  • 
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III.  NUMERICAL  RESULTS 


A.  MONTERE Y-MIAMI  PARABOLIC  EQUATION  GRID  SIZES 

The  discretization  of  the  ocean  environment  is  necessary  and  defined  by  the  mesh 
size  (A r,  A z)  [5].  The  depth  mesh,  A z,  is  used  to  detennine  the  desired  fast  Fourier 
transform  (FFT)  size,  defined  by  nz.  The  relation  between  A z  and  nz  is  defined  as 
nz=2zmax/Az  in  the  SSF  algorithm,  where  zmax  is  the  maximum  computational  depth  and  2 
is  the  factor  to  account  for  the  image  ocean  technique  common  in  SSF  models.  The 
scaling  factor,  rzfact,  is  multiplied  by  the  depth  mesh  to  determine  the  range  step  size,  A r. 
We  get  a  unifonn  grid  size  when  rzfact=  1,  while  the  range  step  size  becomes  greater  than 
the  depth  grid  size  when  rzfact>  1,  or  vice  versa. 

In  this  chapter,  we  first  apply  the  density-smoothing  approach  and  then  the  hybrid 
SSF/FD  method  by  using  different  step  sizes  and  depth  meshes.  In  each  case,  we  compare 
our  computed  transmission  loss  with  a  benchmark  solution  based  on  the  results  computed 
from  the  normal  mode  model,  Couple07,  by  Richard  Evans  [13].  The  comparison 
between  solutions  is  based  on  both  signal  amplitude  and  the  modulation  structure  in 
range  (which  depends  on  the  phase  of  the  solution).  In  many  cases,  the  amplitude  is 
found  to  agree  well,  but  the  phase  structure  shows  disagreement  with  the  benchmark 
solutions. 

After  the  examination  of  these  two  methods,  we  present  a  comparison  between 
the  two  to  show  the  convergence,  sensitivity,  and  other  differences  of  the  density-contrast 
interface.  Finally,  for  completeness,  we  show  the  effect  of  including  the  higher  order 
correction  term. 

For  all  calculations  presented  in  this  thesis,  the  environment  was  based  on  a 
scenario  originally  defined  by  Hursky  [9],  This  consists  of  a  Pekeris  waveguide  of  depth 
300  m  with  a  water  column  sound  speed  of  1500  m/s  overlying  a  semi-infinite  sediment 
half-space  of  sound  speed  1700  m/s.  Except  for  the  initial  test,  in  which  both  the  water 
and  sediment  densities  equal  1.0  g/cm  ,  the  sediment  density  is  defined  as  1.5  g/cm  . 
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Attenuation  in  the  sediment  is  defined  as  0.5  dB/A.  or  0.294  dB/m/kHz,  while  the  water 
column  is  assumed  to  be  lossless.  Pressure  calculations  were  made  for  a  point  source  at 
depth  180  m  transmitting  a  single  frequency  of  100  Hz.  For  comparison,  a  transmission 
loss  trace  for  a  receiver  depth  of  100  m  was  extracted  out  to  20  km  range. 

B.  NO  DENSITY  DISCONTINUITY  CASE 

As  has  been  stated,  the  SSF  algorithm  is  known  to  be  very  accurate  in  the 
prediction  of  transmission  loss  in  the  absence  of  the  density  discontinuities.  In  this  paper, 
we  confirmed  the  consistency  of  the  SSF  algorithm  with  our  benchmark  Couple07 
solution  for  such  an  environment.  Figure  1  shows  the  comparison  of  these  two  solutions 
out  to  20  km  range.  The  agreement  is  found  to  be  excellent. 
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Figure  1.  No  Density  Discontinuity  vs.  Couple07 
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c. 


DENSITY-SMOOTHING  APPROACH 


In  the  presence  of  density  discontinuities  across  boundaries,  the  SSF  algorithm 
utilizes  the  density-smoothing  approach,  as  defined  in  the  previous  section.  The  mixing 
length,  L,  is  a  free  variable  that  can  be  adjusted.  In  this  work,  we  examined  various 
mixing  lengths:  l/4k0,  l/2ko,  1/ko,  3/2ko,  2/ko,  where  ko  is  the  reference  wavenumber.  The 
best  results  were  obtained  when  L  was  set  to  1/ko.  We  also  tried  different  depth  meshes, 
A z,  by  adjusting  nz,  and  range  step  sizes,  A r,  by  changing  the  scaling  factor  rzfact.  We 
acquired  an  optimum  stable,  convergent  solution  with  grid  size  Az=A/15  and  rzfact=  5, 
which  corresponded  to  Az=l  m  with  an  FFT  size  of  «z=4096,  and  range  step  size  Ar=5 
m,  or  A/3.  In  order  to  see  the  effects  of  the  depth  mesh  and  range  step  size  on  the  results 
using  the  density-smoothing  approach,  Figure  2  shows  different  rzfact  values  with  fixed 
Az=A/15,  and  Figure  4  shows  different  A z  values  with  fixed  rzfact= 5.  Each  result  is 
compared  to  our  benchmark  solution. 
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Figure  2. 


Density-Smoothing  Approach,  rzfact=\,  2,  5,  10,  20 


Figure  2  shows  the  transmission  loss  comparisons  with  fixed  value  of  A z  =)J  1 5 . 
The  scaling  factor  rzfact  is  chosen  to  be  1,  2,  5  in  the  first  panel  and  5,  10,  20  in  the 
second  panel.  If  we  examine  the  plot,  the  smoothing  approach  matches  the  benchmark 
solution  very  well  for  the  first  several  km  for  all  rzfact  values  considered.  However, 
between  4  and  10  km,  a  phase  error  appears  and  becomes  significant,  especially  after  10 
km.  It  is  worth  noting  that,  for  Az  <  /i  / 1 0  ■  convergence  appears  to  occur  for  rzfact=  5, 
with  only  very  minor  changes  for  smaller  rzfact  values. 
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Figure  3.  Convergence  for  rzfact=  1,  2,  5,  10,  20 


Figure  3  further  illustrates  the  convergence  of  transmission  loss  results  between 
15  and  20  km.  After  examining  different  rzfact  values,  it  is  obvious  that  we  cannot  see  a 
significant  change  for  rzfact< 5.  However,  when  rzfact>5,  the  transmission  loss  amplitude 
changes  slowly.  For  all  cases,  the  phase  error  remains  consistent  and  is  unaffected  by 
variations  in  rzfact. 
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2. 
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Figure  4.  Density-Smoothing  Approach,  Az=X/2,  a/5,  a/  1 0,  a/1  5,  a/40 


Figure  4  shows  different  values  of  depth  mesh,  A z,  for  a  fixed  value  of  rzfact=  5. 
The  solution  converges  for  A z  <  A/10 ,  and  Az  =  A/15  provides  good  agreement  with  the 
benchmark  solution  for  the  first  several  km.  The  solution  begins  to  degrade  at  shorter 
ranges  for  Az  >  A/5 ,  and  when  Az<a/2,  the  agreement  is  poor  at  all  ranges.  Even  for  the 
best  solution,  however,  the  phase  error  begins  to  appear  at  ranges  beyond  about  8  km. 
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Figure  5.  Convergence  for  Az=)J2,  1/5,  A/ 10,  a/1  5,  a/40 

The  expanded  view  in  Figure  5  shows  a  stable  behavior  for  Az<)J 5 .  For  these 
values,  the  amplitude  stays  constant,  and  we  observe  slight  phase  changes  at  ranges 
beyond  15  km.  However,  smaller  values  of  Az  do  not  improve  the  phase  match  with  the 
benchmark  solution.  In  fact,  for  values  of  Az<)J  1 5 ,  numerical  sensitivities  begin  to 
degrade  the  solution  in  the  form  of  small-scale  fluctuations  in  the  result. 

D.  HYBRID  SPLIT-STEP/FINITE  DIFFERENCE  ALGORITHM 

For  the  hybrid  SSF/FD  approach,  as  with  in  the  density-smoothing  approach,  we 
applied  different  range  step  sizes  and  depth  meshes  to  find  the  optimum  solution.  Using 
the  hybrid  approach,  we  found  improved  agreement  with  our  benchmark  solution  at 
longer  ranges.  We  acquired  the  best  result  by  applying  a  depth  mesh  Az=l/\5=\  m, 
corresponding  to  an  FFT  size  of  nz=4096,  and  rzfact=  10  or  A/^10  m. 
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Figure  6.  Hybrid  Split-Step  Method,  rzfacP=  1,  2,  5,  10,  20 


The  first  panel  in  Figure  6  shows  the  results  for  rzfact=  1,  2,  5,  and  the  second 
panel  shows  them  for  rzfact=  5,  10,  20  with  a  fixed  value  of  Az=)J  1 5 .  It  can  be  concluded 
from  the  plot  that  when  rzfact< 5,  then  the  transmission  loss  amplitude  and  phase  error 
remain  the  same.  The  phase  error  begins  to  appear  after  8  km  but  remains  small  for  all 
rzfact  values. 
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Hybrid  Split-Step/Finite  Difference  Method  Az=A,/15 
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Figure  7.  Convergence  for  rzfact=\,  2,  5,  10,  20 


From  the  expanded  view  in  Figure  7,  an  interesting  observation  is  noted.  First,  as 
in  the  smoothing  approach,  the  solution  appears  to  converge  to  a  stable  result  for 
rzfact< 5.  However,  in  this  case,  a  better  match  with  the  benchmark  result  is  found  just 
prior  to  convergence  when  rzfact=  10.  Larger  values  of  rzfact  give  much  worse  results,  so 
the  improvement  for  rzfact=  10  is  likely  just  a  coincidence.  Therefore,  the  following 
results  are  computed  for  rzfact=  5. 
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2. 
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Figure  8.  Flybrid  Split-Step  Method,  Az=a/5,  a/ 10,  a/1  5,  a/40,  a/60 


Figure  8  shows  the  various  solutions  at  a  fixed  value  of  rzfact=  5  for  A z=  )J  1 5  to 
A z=  a/60.  It  is  obvious  from  the  figure  that  decreasing  the  A z  values  does  not  improve  the 
result  for  rzfact=5  but  introduces  significant  variability  in  the  transmission  loss 
amplitude.  The  results  show  the  best  match  with  Az=A,/15,  where  agreement  with  the 
benchmark  solution  is  good  beyond  10  km.  For  Az  >  A  /  5  (results  computed  but  not 
shown),  agreement  degrades  after  a  few  kilometers,  with  very  poor  agreement  at  all 
ranges  for  Az  «  A, . 
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Figure  9.  Convergence  for  Az=A/5,  A/10,  A/15,  AMO,  A/60 


Figure  9  displays  an  expanded  view  of  the  last  5  km  of  the  simulation  as  the  A z 
values  decreased.  The  A z  values  larger  than  A/ 15  cause  high  phase  errors  and  amplitude 
fluctuations,  while  we  see  good  agreement  with  Az=A/15.  It  seems  that  the  hybrid  method 
is  very  sensitive  to  A z  values  and  show  poor  convergence  in  A z. 

E.  COMPARISON  OF  DENSITY-SMOOTHING  APPROACH  AND  HYBRID 
SPLIT  STEP  ALGORITHM 

The  comparison  of  the  density-smoothing  and  hybrid  split  approaches  gives 
unique  results.  Figure  10  shows  the  optimal  solution  of  the  two  methods  compared  to 
benchmark  solution  Couple07.  Both  methods  give  the  best  result  when  Az=A/15  and 
converge  when  rzfact=  5.  The  agreement  with  the  benchmark  is  consistent  out  to  about  8 
km  for  the  smoothing  approach,  while  this  number  is  almost  18  km  for  the  hybrid 
method. 


25 


Both  the  smoothing  approach  and  the  hybrid  method  tend  to  converge  when  rzfact 
is  around  5.  For  both  methods,  the  results  are  stable  when  rzfact< 5.  In  either  case,  it  is 
recommended  that  rzfact  should  be  between  2  and  5  in  order  to  get  accurate  results. 
Figure  1 1  illustrates  the  last  5  km  of  the  solution  where  hybrid  method  shows  improved 
agreement  with  the  benchmark  solution.  The  transmission  loss  amplitudes  of  these  two 
methods  are  almost  the  same,  but  the  phase  error  of  the  smoothing  approach  is 
significantly  higher  than  the  hybrid  method  at  these  ranges. 

One  of  the  main  differences  between  these  two  methods  is  their  sensitivity  to  A z, 
or  nz.  It  can  be  observed  from  Figure  5  and  Figure  9  that  the  hybrid  method  is  more 
sensitive  than  the  smoothing  approach  since  it  changes  significantly  by  varying  the  Az 
values  and  does  not  appear  to  reach  a  stable,  convergent  solution.  So  while  the  hybrid 
method  appears  to  produce  more  accurate  results,  stability  of  this  method  appears  to 
remain  an  issue. 
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Figure  10.  Smoothing  and  Flybrid  Methods  vs.  Couple07  with  Az=)J  1 5  and 

Ar=5*Az=m 
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Comparison  of  Smoothing  and  Hybrid  Methods 


Figure  1 1 .  Convergence  of  Smoothing  and  Hybrid  Methods  with  Az=kH5  and 

Ar=5*Az=X/3 


F.  COMPARISON  OF  HYBRID  AND  HYBRID  HIGHER  ORDER 
CORRECTION 

Equations  (15)  and  (36)-(43)  define  expressions  that  provide  a  higher  order 
correction  to  the  square  root  operator.  In  Figures  12  and  13,  we  see  the  comparison  of  the 
basic  hybrid  method  and  the  hybrid  higher  order  correction  with  the  benchmark  solution. 
For  our  case,  we  can  conclude  that  higher  order  correction  does  not  improve  the  solution 
any  further. 
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TL(dB  re  1m)  for  Pressure 


Comparison  of  Hybrid  and  Hybrid  Higher  Order  Correction 


Figure  12.  Comparison  of  Hybrid  and  Hybrid  Higher  Order  Correction  with  Az=)J  1 5 

and  A r=  1  *  Az=)J  1 5 
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TL(dB  re  1m)  for  Pressure 


Comparison  of  Hybrid  and  Hybrid  Higher  Order  Correction 


Figure  1 3 .  Comparison  of  Hybrid  and  Hybrid  Higher  Order  Correction  with  Az=A,/ 1 5 

and  A r=  1  *Az=)J  1 5 
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IV.  SUMMARY 


The  MMPE  model,  based  on  the  SSF  marching  algorithm,  has  been  used  to 
predict  the  sound  propagation  in  deep  and  shallow  water  environments  for  many  years. 
Many  studies  have  confirmed  MMPE’s  accuracy  in  a  wide  variety  of  environments. 
However,  in  the  presence  of  boundary  density  discontinuities,  phase  errors  have  been 
known  to  appear  due  to  bottom  interactions.  This  degrades  the  model’s  long-range 
accuracy.  In  this  thesis,  we  compared  results  of  the  split-step  density-smoothing  and 
hybrid  split-step/finite  difference  methods  to  treat  the  density  discontinuity  at  the  water- 
bottom  interface.  Of  particular  interest  was  the  potential  reduction  in  the  phase  errors 
using  the  hybrid  approach,  which  could  then  be  easily  implemented  into  existing  versions 
for  the  MMPE  with  the  goal  of  improving  the  model’s  long  range  accuracy. 

For  the  simple  environment  examined  in  this  thesis,  the  SSF  was  found  to  give 
results  consistent  with  a  benchmark  solution  out  to  several  kilometers,  but  beyond  that, 
the  phase  errors  accumulated  and  became  significant  after  10  km.  In  our  examination,  we 
obtained  the  optimum  result  by  defining  a  depth  mesh,  Az=k/15.  We  also  found  that 
decreasing  the  depth  mesh  does  not  decrease  the  phase  error  significantly,  and  that  the 
method  tends  to  converge  when  A z  is  between  )J  1 0  and  A/15.  Similarly,  the  range  step 
size  plays  an  important  role  in  determining  the  best  solution.  The  optimum  result  for  that 
approach  is  obtained  by  a  range  step  size  Ar=AV3,  corresponding  to  scaling  factor, 
rzfact= 5. 

The  hybrid  split-step/finite  difference  method  resulted  in  improved  solutions 
compared  to  our  benchmark.  This  approach  produced  results  that  closely  matched  the 
benchmark  results  computed  using  Couple07.  The  best  result  was  obtained  with  a  depth 
mesh  Az=A715  and  range  step  size  Ar=2)J3  corresponding  to  rzfact=  10.  However, 
convergence  was  achieved  at  rzfact=  5.  It  is  worth  noting  that  both  methods  produced 
optimal  results  for  the  same  mesh  sizes. 

For  both  methods,  we  found  that  for  a  fixed  depth  mesh  of  Az « Z/1 5 ,  the 
solution  converged  to  a  reasonably  stable  state  for  A r  <  5Az  =  Z/3  .  However,  for  a  fixed 
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range  step  size  of  A  r  =  A/3  ,  only  the  split-step  smoothing  method  appeared  to  converge 
to  a  stable  state  for  decreasing  values  of  Az.  The  hybrid  method  appeared  to  reach  an 
optimal  solution  at  roughly  the  same  depth  mesh  scale  as  the  split-step  smoothing 
method,  Az  «  A/15  ,  but  appeared  to  degrade  for  much  smaller  values  of  Az.  This  suggests 
that  the  finite  difference  approach  to  the  density  discontinuity  is  much  more  sensitive  to 
the  choice  of  Az.  While  the  solution  is  improved,  this  lack  of  stability  should  be  treated 
with  care.  Future  work  should  investigate  the  differences  in  sensitivity  and  convergence 
between  these  two  methods  in  more  detail. 

Finally,  we  evaluated  the  influence  of  the  higher  order  correction  in  the  hybrid 
method.  For  this  simple  environment,  it  did  not  improve  the  solution  significantly.  Due  to 
the  added  complexity  of  this  correction  and  the  lack  of  significant  improvement,  its 
inclusion  in  future  updates  of  the  MMPE  model  is  not  recommended. 
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APPENDIX.  MATLAB  CODES  FOR  HYBRID  METHOD 


This  Matlab  program  performs  a  simple,  range-independent  calculation  of  a 
Pekeris  waveguide  problem  as  defined  by  Paul  Hursky.  The  calculation  is  based  on  a 
split-step  Fourier  method  which,  when  using  the  density-smoothing  approach  defined  by 
Tappert  and  utilized  in  MMPE,  shows  phase  errors  that  accumulate  with  range.  In  this 
code,  the  hybrid  method  of  Yevick  and  Thomson  is  implemented  in  order  to  test  the 
correction  of  the  phase  error  by  eliminating  the  smoothing  approach. 

%%  PARAMETERS 
format  long 

zbot=300;  %  bottom  interface  depth 

D=4*zbot;  %  max  computational  depth  of  real  ocean 

zs=180;  %  source  depth 

freq=100;  %  CW  signal 

cw=1500;  %  sound  speed  in  the  ocean 

cb=1700;  %  sound  speed  in  the  bottom 

c0=1500;  %  reference  sound  speed 

lambda=cO/freq;  %  reference  wavelength 

dz=lambda/10;  %  first  estimate  of  dz 

NZ=round(2  *  D/dz) ;  %  first  estimate  of  NZ 

NZ=nextpow2(NZ);  %  round  up  to  next  power  of  2  for  FFT  size 

NZ=2.ANZ; 

D=(NZ/2)*dz;  %  max  computational  depth 
zl=(l/2)*D;  %  depth  to  begin  spatial  domain  filtering 
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nb=floor(zbot/dz)+l;  %  depth  index  of  bottom  interface 
rzfact=2;  %  scaling  factor  between  dr  and  dz  ( 
dr=rzfact*dz;  %  range  step  size 

dkz=pi/D;  %  vertical  wavenumber  sampling 

kO=(2*pi*freq)/cO;  %  reference  wavenumber 

kl=(3/4)*k0;  %  wavenumber  to  begin  wavenumber  domain  filtering 

R0=  1 ;  %  reference  range 

Rmax=20000;  %  max  computational  range 

NR=floor(Rmax/dr);  %  number  of  range  steps  in  calculation 

%  Density  Parameters 

rho_w=1000;  %  density  in  water 

rho_b=1500;  %  density  in  bottom 

%  %Preallocate  arrays/matrices 

kz=zeros(l,NZ); 

x=zeros(  1  ,NZ); 

A=zeros(l,NZ); 

B=zeros(l,NZ); 

psi=zeros(l,NZ); 

FK=ones(l,NZ); 

Top=zeros(l,NZ); 

PROPK=zeros(  1  ,NZ); 
z=zeros(l,NZ/2); 

FZ=ones(  1  ,NZ/2); 
cz=zeros(l,NZ/2); 
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rho=zeros(  1  ,NZ/2); 
rho_m=zeros(  1  ,NZ/2); 
rho_t=zeros(  1  ,NZ/2); 
rho_p=zeros(  1  ,NZ/2); 
n_prime_sq=zeros(  1  ,NZ/2); 
atten=zeros(  1  ,NZ/2); 
atten_e=zeros(  1  ,NZ/2); 

U  op=zeros(  1  ,NZ/2); 
PROPZ=zeros(l,NZ/2); 
press=zeros(NZ,NR); 
r=zeros(l,NR); 

TL=zeros(NZ,NR); 

%  Hybrid  Method  Parameters 

se=(kO*dz).A2; 

delta=kO*dr; 

alpha=(  1  /4)  *  ( 1  - 1  i  *  delta) ; 

beta=(l/4)*(l+li*delta); 

q  1 =(rho_w-rho_b  )/(3  *  rho_w+rho_b ) ; 

q2=(rho_b-rho_w)/(3*rho_b+rho_w); 

q3 =(rho_w-rho_b )/ (rho_w+3  *  rho_b ) ; 

q4=(rho_b-rho_w)/(rho_b+3*rho_w); 

al=(l-(alpha/se)*ql); 

a2=(alpha/se)*ql; 

a3=(alpha/se)*q4; 
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a4=(l-(alpha/se)*(q3+q4)); 

a5=(alpha/se)*q3; 

a6=(alpha/se)*q2; 

a7=(l-(alpha/se)*q2); 

bl=(l-(beta/se)*ql); 

b2=(beta/se)*ql; 

b3=(beta/se)*q4; 

b4=(  1  -(beta/se)*(q3+q4)); 

b5=(beta/se)*q3; 

b6=(beta/se)*q2; 

b7=(l-(beta/se)*q2); 

%  factors  used  in  algebraic  solution  of  tridiagonal  matrix 
factO=  1  ,/(a4-a2  *  a3/a  1  -a5  *  a6/ a7) ; 
factl=(b3-bl*a3/al)*fact0; 
fact2=(b4-b2  *  a3/a  1  -b6  *a5/a7)  *factO ; 
fact3=(b5-b7*a5/a7)*fact0; 

%  produce  r=0  in  kz-domain 

kz=[-NZ/2:(NZ/2-l)]*dkz; 

x=kz/kO; 

A=(l-x.A2).A(-l/4); 

B=-li*NZ*dkz*sqrt(2/(pi*kO))*sin(kz*zs); 

psi=A.*B; 
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%  Wavenumber  domain  filter 


for  ii=l:NZ; 

if  -kO<kz(ii)  &&  kz(ii)<-kl  ; 

FK(ii)=(cos((2  *pi/kO)  *  (kz(ii)-k  1 ))) .  A2 ; 
elseif  kl<kz(ii)  &&  kz(ii)<kO; 

FK(ii)=(cos((2  *pi/kO)  *  (kz(ii)-k  1 ))) .  A2 ; 
elseif  kz(ii)<-kO; 

FK(ii)=0; 
elseif  kO<kz(ii); 

FK(ii)=0; 

end 

end 

FK=fftshift(FK); 

psi=fftshift(psi).*FK;  %  starting  field  with  wavenumber  filter  applied  in  FFT 

order 

psi=ifft(psi);  %  transform  to  z-domain 
%  Diffraction  propagator  (kz-domain) 

T  op=  1  -sqrt(  1  -x.  A2); 

PROPK=exp(-li*kO*dr*Top); 

PROPK=fftshift(PROPK);  %  put  in  FFT  order  for  application  in  SSF  algorithm 
%  Depth 

z=[0:NZ/2-l]*dz; 
for  IZ=l:NZ/2; 
if  0<=z(IZ)  &&  z(IZ)<=zl; 
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FZ(IZ)=1; 


elseif  zl<z(IZ)  &&  z(IZ)<=(D-dz); 

FZ(IZ)=(cos((pi/ (2  *  (D-z  1 )))  *  (z(IZ)-z  1 ))).  A2; 
end 
end 

%  Sound  Speed  Profile  (real  ocean  part  of  z-domain) 

cz=size(z); 

cz(l:nb-l)=cw; 

cz(nb)=0 . 5  *  (cw+cb) ; 

cz(nb+l  :NZ/2)=cb; 

n_prime_sq=(cO./cz).A2;  %  For  density  smoothing  approach  use  effective  index  of 
refraction  instead 

%  Density  Profile  (real  ocean  part  of  z-domain) 

rho=size(z); 

rho(  1  :nb- 1  )=rho_w; 

rho(nb)=0 . 5  *  (rho_w+rho_b) ; 

rho(nb+l  :NZ/2)=rho_b; 

%  Bottom  attenuation 

atten(nb)=0.5*(0.5);  %  in  dB/lambda  -  applying  Heaviside  step  function 
intermediate  value 

atten(nb+l:NZ/2)=0.5;  %  in  dB/lambda 
atten_e=atten.*(freq./cz)/8.686;  %  in  1/m 
%  Lens  propagator  (real  ocean  part  of  z-domain) 

Uop=(  1  -sqrt(n_prime_sq)); 
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PROPZ=exp(-li*kO*dr*Uop).*exp(-atten_e*dr); 

%  Range  increament 
for  IR=1:NR; 

psi(l:NZ/2)=PROPZ.*psi(l:NZ/2);  %  lens  propagator  to  real  ocean 
psi(  1  :NZ/2)=FZ.  *psi(  1  :NZ/2);  %  z-domain  filter  to  real  ocean 
psi(NZ/2+2:NZ)=-psi(NZ/2:-l:2);  psi(NZ/2+l)=0;  %  odd  symmetry 
psi=fft(psi);  %  transfonn  to  kz-domain  -  output  in  FFT  order 
psi=psi.*PROPK;  %  diffraction  propagator 
psi=psi.*FK;  %  kz-domain  filter 

psi=ifft(psi);  %  transform  to  z-domain  -  output  in  FFT  order 

%  Flybrid  method  -  apply  in  real  ocean  only  (odd  symmetry  applied  after  lens 
propagator) 

if  IR==  1 

disp('Using  hybrid  finite  difference  method  for  density') 
end 

tmp=fact  1  *psi(nb- 1  )+fact2  *psi(nb)+fact3  *psi(nb+ 1 ); 
psi(nb-l)=b  1/al  *psi(nb-  l)+b2/al  *psi(nb)-a2/al  *tmp; 
psi(nb+ 1  )=b6/ a7  *psi(nb)+b7/a7  *psi(nb+ 1  )-a6/ a7  *  tmp ; 
psi(nb)=tmp; 

%  Compute  Pressure  Field 
r(lR)=lR*dr; 

p res s( : ,  I R  )=p s  i  *  sq rtf  R0/r(  I R ) ) ;  %  store  pressure  field  in  FFT  order 
end 
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%  Transmission  Loss 


TL=-20*log  1 0(abs(press(  1  :NZ/2, :))); 
zout=100; 

nout=find(z>=zout,  1  ); 
if  (zout-z(nout- 1  ))<(z(nout)-zout) 
nout=nout-l; 
end 
figure 

plot(r/ 1 000,TL(nout, :),'r');axis  ij ; 
xlim  ( [0  20]) 
ylim  ([30  110]) 
hold  on 

plot(rng,tlz,’k’);axis  ij; 
xlim  ([0  20]) 
ylim  ([30  110]) 

legend(’PE-SSF-YT-lstOrder’,'Couple07') 

title(['dr-  num2str(dr)  ’m,  lambda/dz-  num2str(lambda/dz)  c0=  num2str(c0) 
'm/s,  rhob=  num2str(rho_b)  'kg/mA3  ]) 
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